library(metafor)
library(plyr)
library(tidyverse)
library(magrittr)
library(foreign)
library(estimatr)
library(emmeans)
library(broom)
library(countrycode)


t1 <- "https://github.com/thomasjwood/kz_multicountry/raw/main/t2.rds" %>% 
  url %>% 
  gzcon %>% 
  readRDS %>% 
  group_by(
    geo_country_code, issue, wave
  ) %>% 
  nest %>% 
  mutate(
    mods = data %>% 
      map(
        \(i)
        
        lm_robust(
          out_num ~ treatment,
          data = i
        )
      ),
    emm = mods %>% 
      map(
        ~emmeans(
          ., 
          consec ~ treatment, 
          adjust = "none"
        )
      )
  )


# estimate meta analytic models -------------------------------------------

t4_2 <- t1 %>%
  filter(wave == "wave 1") %>% 
  ungroup %>% 
  select(geo_country_code, issue, emm) %>% 
  pmap_dfr(
    function(
    geo_country_code, issue, emm
    )
      
      # emm <- t4$emm[[1]]
      
      emm$contrasts %>% 
      tidy %>% 
      mutate(
        geo_country_code, issue, 
        wave = "wave 1",
        geo_country_code = "meta"
      )
  ) %>% 
  nest_by(
    issue, contrast
  ) %>% 
  ungroup %>% 
  mutate(
    rma = data %>%
      map(
        \(i)
        rma(yi = i$estimate, sei = i$std.error)
        )
    )

# overall meta analytics

t4_3 <- t1 %>%
  filter(wave == "wave 1") %>% 
  ungroup %>% 
  select(geo_country_code, issue, emm) %>% 
  pmap_dfr(
    function(
    geo_country_code, issue, emm
    )
      
      emm$contrasts %>% 
      tidy %>% 
      mutate(
        geo_country_code, issue, 
        wave = "wave 1",
        geo_country_code = "meta"
      )
  ) %>% 
  nest_by(
    contrast
  ) %>% 
  ungroup %>% 
  mutate(
    rma = data %>%
      map(
        \(i)
        rma(
          yi = i$estimate, 
          sei = i$std.error
          )
        )
    )

# extracting estimates ----------------------------------------------------

t5 <- t1 %>% 
  ungroup %>% 
  select(
    geo_country_code, wave, issue, emm
  ) %>% 
  pmap_dfr(
    \(geo_country_code, wave, issue, emm)
    emm %>% 
      map_dfr(
        ~tidy(.) %>%
          rename_with(
            ~str_replace(., "contrast", "treatment")
          ),
        .id = "type"  
      ) %>% 
      select(-term, -null.value) %>% 
      mutate(
        geo_country_code, wave, issue
      )
  ) %>% 
  mutate(
    cvar = case_when(
      type %>% 
        equals("emmeans") ~ "Conditional means",
      type %>% 
        equals("emmeans") %>% 
        not  &
        treatment %>% 
        str_detect("Misinformation - ") ~ "Misinformation effect",
      type %>% 
        equals("emmeans") %>% 
        not  &
        treatment %>% 
        str_detect("Misinformation &") ~ "Correction effect"
    ) %>% 
      factor(
        c("Conditional means",
          "Misinformation effect",
          "Correction effect")
      )
  ) %>% 
  filter(
    type == "contrasts"
  ) %>% 
  bind_rows(
    # joining meta analytic estimates
    t4_2 %>% 
      ungroup %>% 
      pmap_dfr(
        function(issue, contrast, data, rma)
          
          rma %>% 
          tidy %>% 
          mutate(
            issue, contrast, 
            wave = data$wave[[1]],
            geo_country_code = data$geo_country_code[[1]]
          )
      ) %>% 
      mutate(
        cvar = case_when(
          contrast %>%
            str_detect("Misinformation &") ~ "Correction effect",
          TRUE ~ "Misinformation effect"
        )  %>% 
          factor(
            c("Conditional means",
              "Misinformation effect",
              "Correction effect")
          ),
        type = "meta"
      ) %>%
      rename(
        treatment = contrast
      ) %>% 
      select(-term)
  ) %>% 
  # joining overall meta
  bind_rows(
    t4_3 %>% 
      ungroup %>% 
      pmap_dfr(
        function(contrast, data, rma)
          
          rma %>% 
          tidy %>% 
          mutate(
            contrast, 
            issue = "",
            wave = data$wave[[1]],
            geo_country_code = "over_meta"
          )
      ) %>% 
      mutate(
        cvar = case_when(
          contrast %>%
            str_detect("Misinformation &") ~ "Correction effect",
          TRUE ~ "Misinformation effect"
        )  %>% 
          factor(
            c("Conditional means",
              "Misinformation effect",
              "Correction effect")
          ),
        type = "over_meta"
      ) %>%
      rename(
        treatment = contrast
      ) %>% 
      select(-term)
  ) %>%
  mutate(
    yl = issue %>% 
      str_detect(
        "country "
      ) %>% 
      ifelse(
        "cs_",
        issue %>% 
          str_sub(-1) %>% 
          str_c("gl_", .)
      ) %>% 
      str_c(geo_country_code)
  ) %>% 
  filter(wave == "wave 1")

t5$labs <- t5$estimate %>%
  round(2) %>% 
  abs %>% 
  equals(0) %>% 
  ifelse(
    t5$estimate %>% 
      round(3) %>%
      str_replace("0\\.", "."),
    t5$estimate %>% 
      round(2) %>%
      str_replace("0\\.", ".")
  )

t5$labs %<>% 
  str_detect("-") %>%
  ifelse(
    t5$labs %>% 
      str_pad(width = 4, side = "right", pad = "0"),
    t5$labs %>% 
      str_pad(width = 3, side = "right", pad = "0")
  )

t5$labs[t5$labs == "000"] <- t5$estimate[t5$labs == "000"] %>% 
  round(6) %>%
  format(scientific = F) %>% 
  str_replace("0\\.", ".")

t5$labs %<>% 
  str_c(
    t5$p.value %>% 
      gtools::stars.pval() %>% 
      str_remove("\\.")
  )


t5$low <- t5$estimate %>% 
  subtract(
    t5$std.error %>% 
      multiply_by(1.96)
  )

t5$hi <- t5$estimate %>% 
  add(
    t5$std.error %>% 
      multiply_by(1.96)
  )



# generate yaxis lab ------------------------------------------------------

t5$yl %<>% 
  factor(
    t5 %>%
      filter(
        wave == "wave 1" &
          cvar == "Correction effect"
      ) %>% 
      arrange(
        issue, cvar, type, issue, desc(estimate)
      ) %>% 
      use_series(yl) %>% 
      str_remove_all("cs_meta|gl_1meta|gl_2meta|gl_over_meta") %>% 
      stringi::stri_omit_empty() %>% 
      c("cs_meta","gl_1meta", "gl_2meta", "gl_over_meta", .) %>% 
      mapvalues("gl_over_meta", "gl_")
  )


t5$issue %<>% 
  mapvalues(
    c("global 1",
      "global 2",
      "country specific",
      ""),
    c("The Gates Foundation is going to make huge profits from COVID-19 vaccines.",
      "mRNA vaccines like Moderna permanently alter DNA, according to Moderna executives.",
      "Country specific misinformation",
      "Overall") %>% 
      str_wrap(17)
  )


t5$issue %<>%
  factor(
    t5 %>% 
      filter(
        cvar %>% 
          str_detect("Correction")
      ) %>% 
      group_by(issue) %>% 
      summarize(mu = estimate %>% mean) %>%
      arrange(mu) %>% 
      use_series(issue)
  )

t5$issue %<>% 
  fct_relevel("Overall", after = 3) 

# plotting ----------------------------------------------------------------

ggplot() +
  geom_vline(
    xintercept = 0,
    linetype = "dotted",
    size = .25
  ) +
  geom_blank(
    aes(
      estimate, yl
    ),
    data = t5
  ) +
  geom_label(
    aes(
      estimate, yl, label = labs,
    ),
    label.size = 0,
    data = t5 %>% 
      filter(
        yl %>% 
          str_detect("meta", negate = T)
      ),
    size = 2.25,
    fill = "grey95",
    position = position_nudge(0, .4)
    ) +
  geom_point(
    aes(
      estimate, yl
    ),
    data = t5 %>% 
      filter(
        yl %>% 
          str_detect("meta", negate = T)
      ),
    fill = "white"
  ) +
  geom_linerange(
    aes(
      xmin = low, xmax = hi, y = yl,
    ),
    data = t5,
    size = .25
  ) +
  geom_point(
    aes(
      estimate, yl
    ),
    data = t5 %>% 
      filter(
        yl %>% 
          str_detect("meta", negate = F)
      ),
    size = 8.5,
    shape = 23,
    fill = "white"
  ) +
  geom_point(
    aes(
      estimate, yl
    ),
    data = t5 %>% 
      filter(
        yl %>% 
          is.na
      ),
    size = 8.5,
    shape = 23,
    fill = "white"
  ) +
  geom_label(
    aes(
      estimate, yl,
      label = labs
    ),
    size = 2.25,
    data = t5 %>% 
      filter(
        yl %>% 
          str_detect("meta")
      ),
    label.size = 0,
    fill = "transparent",
    fontface = "bold",
  ) +
  geom_label(
    aes(
      estimate, yl,
      label = labs
    ),
    size = 2.25,
    data = t5 %>% 
      filter(
        yl %>% 
          is.na
      ),
    label.size = 0,
    fill = "transparent",
    fontface = "bold",
  ) +
  facet_grid(
    issue ~ cvar,
    scale = "free",
    space = "free"
    # labeller = label_wrap_gen(width = 20)
  ) +
  scale_y_discrete(
    breaks = t5$yl %>% 
      levels,
    labels = t5$yl %>% 
      levels %>% 
      str_remove_all("gl_1|gl_2|cs_") %>% 
      str_replace_all("meta", "Meta Analysis"),
    expand = expansion(add = c(.9, .8))
  ) +
  scale_x_continuous(
    breaks = seq(-.6, .2, by = .2),
    labels = c("-.6", "-.4", "-.2", "0", ".2")
  ) +
  labs(x = "",
       y = "") 